///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/fqmc/vmc_walker.cpp
///
/// \brief   implementation for a class for proposing moves based on distances between orbitals
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include<string>
#include<vector>
#include<set>
#include<cmath>
#include<algorithm>
#include<numeric>
#include<iostream>

#include<formic/fqmc/distance_proposer.h>
#include<formic/random/random.h>
#include<formic/numeric/numeric.h>
#include<formic/wfn/configuration.h>
#include<formic/wfn/symmetry.h>
#include<formic/input/base.h>
#include<formic/timing/timing.h>
#include<formic/lapack/interface.h>

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief Constructs the proposing object.
///
/// \param[in]      userinp  the user's input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////
formic::fqmc::DistProposer::DistProposer(const formic::InputBase & userinp)
  : _norb(userinp.get<int>("nao")),
    _na(userinp.get<int>("nalpha")),
    _nb(userinp.get<int>("nbeta")),
    _nirr(userinp.get<std::vector<int> >("opirr").size()),
    _sig2( formic::square_norm(userinp.get<double>("dist_move_sd")) ),
    _2e_prob(userinp.get<double>("te_move_prob")),
    _conf(userinp),
    _dpt(),
    _orb_irrep(userinp.get<std::vector<int> >("orb_irrep")),
    _orbc( userinp.get<std::vector<double> >("orb_coords") ),
    _np(0),
    _plist( 12 * userinp.get<int>("nao") * userinp.get<int>("nao"),  -1 ),
    _dw(    6 * userinp.get<int>("nao") * userinp.get<int>("nao"), 0.0 ),
    _n_all_aa_hole_pairs(0),
    _n_all_bb_hole_pairs(0),
    _n_all_ab_hole_pairs(0),
    _all_aa_hole_pairs( 3 * userinp.get<int>("nao") * userinp.get<int>("nao"), -1 ),
    _all_bb_hole_pairs( 3 * userinp.get<int>("nao") * userinp.get<int>("nao"), -1 ),
    _all_ab_hole_pairs( 3 * userinp.get<int>("nao") * userinp.get<int>("nao"), -1 ),
    _aa_hole_pair_irrep_count(userinp.get<std::vector<int> >("opirr").size(), 0),
    _bb_hole_pair_irrep_count(userinp.get<std::vector<int> >("opirr").size(), 0),
    _ab_hole_pair_irrep_count(userinp.get<std::vector<int> >("opirr").size(), 0),
    _n_all_aa_part_pairs(0),
    _n_all_bb_part_pairs(0),
    _n_all_ab_part_pairs(0),
    _all_aa_part_pairs( 3 * userinp.get<int>("nao") * userinp.get<int>("nao"), -1 ),
    _all_bb_part_pairs( 3 * userinp.get<int>("nao") * userinp.get<int>("nao"), -1 ),
    _all_ab_part_pairs( 3 * userinp.get<int>("nao") * userinp.get<int>("nao"), -1 ),
    _aa_part_pair_irrep_count(userinp.get<std::vector<int> >("opirr").size(), 0),
    _bb_part_pair_irrep_count(userinp.get<std::vector<int> >("opirr").size(), 0),
    _ab_part_pair_irrep_count(userinp.get<std::vector<int> >("opirr").size(), 0),
    _chosen_pair_irrep_counts(userinp.get<std::vector<int> >("opirr").size(), 0),
    _dw_mat(userinp.get<int>("nao") * userinp.get<int>("nao"), 0.0)
{

  // get the point group's direct product table
  formic::group_direct_product_table(userinp.get<std::string>("point_group"), _dpt);

  // check number of orbital coordinates
  if ( _orbc.size() != 3*_norb )
    throw formic::Exception("wrong number of orbital coordinates in DistProposer::DistProposer");

  // compute the elements of the distance weighting matrix
  for (int p = 0; p < _norb; p++)
  for (int q = 0; q <= p; q++) {
    const double weight = std::exp( -(   formic::square_norm( _orbc.at(3*p+0) - _orbc.at(3*q+0) )
                                       + formic::square_norm( _orbc.at(3*p+1) - _orbc.at(3*q+1) )
                                       + formic::square_norm( _orbc.at(3*p+2) - _orbc.at(3*q+2) ) ) / _sig2 );
    _dw_mat.at( _norb * p + q ) = weight;
    _dw_mat.at( _norb * q + p ) = weight;
  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  counts how many out of the supplied orbitals are in each irrep
///
/// \param[in]      norb        the number of orbitals
/// \param[in]      nirr        the number of irreducible representations
/// \param[in]      orbirr      the irreducible representation of each orbital
/// \param[in]      orbs        the list of orbitals from which to construct the pairs
/// \param[out]     irrcnt      size == nirr, on exit, for each irrep, the number of orbitals
///                             with that irrep
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::fqmc::DistProposer::count_irreps(const int norb,
                                              const int nirr,
                                              const int * const orbirr,
                                              const int * const orbs,
                                              int * const irrcnt) {

  // zero the irrep counts
  for (int i = 0; i < nirr; i++)
    irrcnt[i] = 0;

  // count the number of orbs in each irrep
  for (int i = 0; i < norb; i++)
    irrcnt[orbirr[orbs[i]]] += 1;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  get a list of all pairs and their irreps that can be made from the supplied orbital
///         list, as well as an array counting how many pairs have each irrep
///
/// \param[in]      norb        the number of orbitals
/// \param[in]      nirr        the number of irreducible representations
/// \param[in]      dpt         the direct product table giving the irreps of orbital products
/// \param[in]      orbirr      the irreducible representation of each orbital
/// \param[in]      orbs        the list of orbitals from which to construct the pairs
/// \param[out]     npair       on exit, the number of pairs in the pair array
/// \param[out]     pairs       size == 3*norb*(norb-1)/2, on exit, the orbital pairs and their
///                             irreps
/// \param[out]     irrcnt      size == nirr, on exit, for each irrep, the number of pairs with
///                             that irrep
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::fqmc::DistProposer::prep_pair_array_and_irrep_counts_aa(const int norb,
                                                                     const int nirr,
                                                                     const int * const dpt,
                                                                     const int * const orbirr,
                                                                     const int * const orbs,
                                                                     int & npair,
                                                                     int * const pairs,
                                                                     int * const irrcnt) {

  // zero the irrep counts
  for (int i = 0; i < nirr; i++)
    irrcnt[i] = 0;

  // get the number of pairs
  npair = norb * ( norb - 1 ) / 2;

  // construct the pair array and irrep counts
  int * pair_ptr = pairs;
  for (int i = 1; i < norb; i++) {
    const int orbi = orbs[i];
    const int * const dpt_row = dpt + orbirr[orbi] * nirr;
    for (int j = 0; j < i; j++) {
      const int orbj = orbs[j];
      const int irr = dpt_row[ orbirr[orbj] ];
      irrcnt[irr] += 1;
      pair_ptr[0] = irr;
      pair_ptr[1] = orbi;
      pair_ptr[2] = orbj;
      pair_ptr += 3;
    }
  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  get a list of all pairs and their irreps that can be made from the supplied
///         alpha and beta orbital lists, as well as an array counting how many pairs
///         have each irrep
///
/// \param[in]      na          the number of alpha orbitals
/// \param[in]      nb          the number of beta  orbitals
/// \param[in]      nirr        the number of irreducible representations
/// \param[in]      dpt         the direct product table giving the irreps of orbital products
/// \param[in]      orbirr      the irreducible representation of each orbital
/// \param[in]      orba        the list of alpha orbitals from which to construct the pairs
/// \param[in]      orbb        the list of beta  orbitals from which to construct the pairs
/// \param[out]     npair       on exit, the number of pairs in the pair array
/// \param[out]     pairs       size == 3*na*nb.  On exit, the orbital pairs and their
///                             irreps
/// \param[out]     irrcnt      size == nirr.  On exit, for each irrep, the number of pairs with
///                             that irrep
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::fqmc::DistProposer::prep_pair_array_and_irrep_counts_ab(const int na,
                                                                     const int nb,
                                                                     const int nirr,
                                                                     const int * const dpt,
                                                                     const int * const orbirr,
                                                                     const int * const orba,
                                                                     const int * const orbb,
                                                                     int & npair,
                                                                     int * const pairs,
                                                                     int * const irrcnt) {

  // zero the irrep counts
  for (int i = 0; i < nirr; i++)
    irrcnt[i] = 0;

  // get the number of pairs
  npair = na * nb;

  // construct the pair array and irrep counts
  int * pair_ptr = pairs;
  for (int i = 0; i < na; i++) {
    const int orbi = orba[i];
    const int * const dpt_row = dpt + orbirr[orbi] * nirr;
    for (int j = 0; j < nb; j++) {
      const int orbj = orbb[j];
      const int irr = dpt_row[ orbirr[orbj] ];
      irrcnt[irr] += 1;
      pair_ptr[0] = irr;
      pair_ptr[1] = orbi;
      pair_ptr[2] = orbj;
      pair_ptr += 3;
    }
  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Determines all the hole-hole pairs and particle-particle pairs and their irreps for
///         the currently stored configuration
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::fqmc::DistProposer::prep_pairs_and_irreps() {

  // determine hole-hole pairs and their irreps
  this->prep_pair_array_and_irrep_counts_aa(_conf.a_unocc_to_site().size(),
                                            _nirr,
                                            &_dpt[0],
                                            &_orb_irrep[0],
                                            &_conf.a_unocc_to_site()[0],
                                            _n_all_aa_hole_pairs,
                                            &_all_aa_hole_pairs[0],
                                            &_aa_hole_pair_irrep_count[0]);
  this->prep_pair_array_and_irrep_counts_aa(_conf.b_unocc_to_site().size(),
                                            _nirr,
                                            &_dpt[0],
                                            &_orb_irrep[0],
                                            &_conf.b_unocc_to_site()[0],
                                            _n_all_bb_hole_pairs,
                                            &_all_bb_hole_pairs[0],
                                            &_bb_hole_pair_irrep_count[0]);
  this->prep_pair_array_and_irrep_counts_ab(_conf.a_unocc_to_site().size(),
                                            _conf.b_unocc_to_site().size(),
                                            _nirr,
                                            &_dpt[0],
                                            &_orb_irrep[0],
                                            &_conf.a_unocc_to_site()[0],
                                            &_conf.b_unocc_to_site()[0],
                                            _n_all_ab_hole_pairs,
                                            &_all_ab_hole_pairs[0],
                                            &_ab_hole_pair_irrep_count[0]);

  // determine particle-particle pairs and their irreps
  this->prep_pair_array_and_irrep_counts_aa(_conf.a_occ_to_site().size(),
                                            _nirr,
                                            &_dpt[0],
                                            &_orb_irrep[0],
                                            &_conf.a_occ_to_site()[0],
                                            _n_all_aa_part_pairs,
                                            &_all_aa_part_pairs[0],
                                            &_aa_part_pair_irrep_count[0]);
  this->prep_pair_array_and_irrep_counts_aa(_conf.b_occ_to_site().size(),
                                            _nirr,
                                            &_dpt[0],
                                            &_orb_irrep[0],
                                            &_conf.b_occ_to_site()[0],
                                            _n_all_bb_part_pairs,
                                            &_all_bb_part_pairs[0],
                                            &_bb_part_pair_irrep_count[0]);
  this->prep_pair_array_and_irrep_counts_ab(_conf.a_occ_to_site().size(),
                                            _conf.b_occ_to_site().size(),
                                            _nirr,
                                            &_dpt[0],
                                            &_orb_irrep[0],
                                            &_conf.a_occ_to_site()[0],
                                            &_conf.b_occ_to_site()[0],
                                            _n_all_ab_part_pairs,
                                            &_all_ab_part_pairs[0],
                                            &_ab_part_pair_irrep_count[0]);

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  extract a list of pairs that have allowed irreps
///
/// \param[in]      s1          the spin of the 1st orbital in a pair
/// \param[in]      s2          the spin of the 2nd orbital in a pair
/// \param[in]      npair       the number of candidate pairs
/// \param[in]      nirr        the number of irreducible representations
/// \param[in]      pairs       size 3*npair.  The list of candidate pairs, with each pair stored
///                             as three ints (irrep, orb1, orb2), where orb1/orb2 are spatial
///                             orbital indices
/// \param[in]      irrcnt      size nirr.  Array giving which irreps are allowed.  An irrep is
///                             allowed if its value in this array is greater than zero.
/// \param[out]     allowed     size 3*npair.  On exit, the allowed pairs, with each pair stored
///                             as three ints (irrep, orb1, orb2), where orb1/orb2 are spin
///                             orbital indices
///
/// \return  The number of allowed pairs
///
///////////////////////////////////////////////////////////////////////////////////////////////////
int formic::fqmc::DistProposer::get_irrep_allowed_pairs(const int s1,
                                                        const int s2,
                                                        const int npair,
                                                        const int nirr,
                                                        const int * const pairs,
                                                        const int * const irrcnt,
                                                        int * const allowed) {

  int p = 0;
  const int three_times_npair = 3*npair;
  for (int i = 0; i < three_times_npair; i += 3) {
    if ( irrcnt[ pairs[i+0] ] > 0 ) {
      allowed[p+0] = pairs[i+0];
      allowed[p+1] = 2*pairs[i+1]+s1;
      allowed[p+2] = 2*pairs[i+2]+s2;
      p += 3;
    }
  }
  return p / 3;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  adds pairs of sites satisfying the given conditions to the list of pairs
///
/// \param[in]      s1       the 1st sites's spin
/// \param[in]      s2       the 2nd sites's spin
/// \param[in]      occ      the required site occupation
/// \param[in]      allwi    the allowed irreducible representations of the pairs' direct products
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::fqmc::DistProposer::get_pairs(const int s1,
                                           const int s2,
                                           const int occ,
                                           const std::set<int> & allwi) {

  if ( s1 > s2 )
    throw formic::Exception("DistProposer::get_pairs expects alpha before beta");

  // for same spin we loop over distinct spatial site pairs,
  // while for opposite spin we loop over all spatial site pairs
  for (int i = 0; i < _norb; i++) {
    for (int j = ( s1 == s2 ? i + 1 : 0 ); j < _norb; j++) {
      const int x1 = 2*i+s1;
      const int x2 = 2*j+s2;
      if ( _conf[x1] == occ && _conf[x2] == occ ) {
        if ( allwi.count( this->orb_pair_dp(i,j) ) ) {
          _plist[2*_np+0] = x1;
          _plist[2*_np+1] = x2;
          _np++;
        }
      }
    }
  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  computes a weight for each orbital in _plist based on the distance from the
///         given orbital
///
/// \param[in]      o1       the given spin orbital
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::fqmc::DistProposer::fast_get_distance_weights(const int o1) {

  // get un-normalized weights
  const double * const dw_row = &_dw_mat.at(0) + _norb * (o1/2);
  for (int p = 0; p < _np; p++)
    _dw[p] = dw_row[_plist[p]/2];

  // get sum of un-normalized weights
  const double weight_sum = std::accumulate(_dw.begin(), _dw.begin()+_np, 0.0);

  // normalize the weights
  formic::xscal(_np, 1.0 / weight_sum, &_dw[0], 1);

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  computes a weight for each pair in _plist based on the distances between the pair's
///         orbitals with the given orbitals
///
/// \param[in]      o1       the 1st spin orbital
/// \param[in]      o2       the 2nd spin orbital
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::fqmc::DistProposer::fast_get_distance_weights(const int o1, const int o2) {

  // get un-normalized weights
  if ( o1 % 2 == o2 % 2 ) {

    const double * const dw_row1 = &_dw_mat.at(0) + _norb * (o1/2);
    const double * const dw_row2 = &_dw_mat.at(0) + _norb * (o2/2);
    const int * pair_ptr = &_plist[0];
    for (int p = 0; p < _np; p++) {
      const int x1 = pair_ptr[1]/2;
      const int x2 = pair_ptr[2]/2;
      const double w1 = dw_row1[x1] * dw_row2[x2];
      const double w2 = dw_row1[x2] * dw_row2[x1];
      //const bool w1_gt_w2 = ( w1 > w2 );
      //_dw[p] = int( w1_gt_w2 ) * w1 + int( !w1_gt_w2 ) * w2;
      _dw[p] = std::max(w1, w2);
      pair_ptr += 3;
    }

  } else {

    const double * const dw_row1 = &_dw_mat.at(0) + _norb * (o1/2);
    const double * const dw_row2 = &_dw_mat.at(0) + _norb * (o2/2);
    const int * pair_ptr = &_plist[0];
    for (int p = 0; p < _np; p++) {
      const int x1 = pair_ptr[1]/2;
      const int x2 = pair_ptr[2]/2;
      _dw[p] = dw_row1[x1] * dw_row2[x2];
      pair_ptr += 3;
    }

  }

  // get sum of un-normalized weights
  const double weight_sum = std::accumulate(_dw.begin(), _dw.begin()+_np, 0.0);

  // normalize the weights
  formic::xscal(_np, 1.0 / weight_sum, &_dw[0], 1);

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  computes a weight for each site in the stored list based on the distance from the
///         given orbital
///
/// \param[in]      o1       the orbital from which distances will be measured
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::fqmc::DistProposer::get_distance_weights(const int o1) {

  // get un-normalized weights
  for (int p = 0; p < _np; p++) {
    const int x1 = _plist[p]/2;
    _dw[p] = std::exp( -(   formic::square_norm( _orbc[3*o1+0] - _orbc[3*x1+0] )
                          + formic::square_norm( _orbc[3*o1+1] - _orbc[3*x1+1] )
                          + formic::square_norm( _orbc[3*o1+2] - _orbc[3*x1+2] ) ) / _sig2 );
  }

  // get sum of un-normalized weights
  const double weight_sum = std::accumulate(_dw.begin(), _dw.begin()+_np, 0.0);

  // normalize the weights
  for (int p = 0; p < _np; p++)
    _dw[p] /= weight_sum;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  computes a weight for each pair in the pair list based on the distance from the
///         given orbitals
///
/// \param[in]      o1       the 1st orbital
/// \param[in]      o2       the 2nd orbital
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::fqmc::DistProposer::get_distance_weights(const int o1, const int o2) {

  // get un-normalized weights
  for (int p = 0; p < _np; p++) {
    const int x1 = _plist[2*p+0]/2;
    const int x2 = _plist[2*p+1]/2;
    _dw[p] = std::exp( -(   formic::square_norm( _orbc[3*o1+0] - _orbc[3*x1+0] )
                          + formic::square_norm( _orbc[3*o1+1] - _orbc[3*x1+1] )
                          + formic::square_norm( _orbc[3*o1+2] - _orbc[3*x1+2] )
                          + formic::square_norm( _orbc[3*o2+0] - _orbc[3*x2+0] )
                          + formic::square_norm( _orbc[3*o2+1] - _orbc[3*x2+1] )
                          + formic::square_norm( _orbc[3*o2+2] - _orbc[3*x2+2] )
                        ) / _sig2 );
  }

  // get sum of un-normalized weights
  const double weight_sum = std::accumulate(_dw.begin(), _dw.begin()+_np, 0.0);

  // normalize the weights
  for (int p = 0; p < _np; p++)
    _dw[p] /= weight_sum;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  stores all particles with valid matching holes in the particle list
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::fqmc::DistProposer::get_valid_particles(const int b) {

//  if ( b < 0 ) {
//
//    _np = 0;
//    this->get_valid_particles(0);
//    this->get_valid_particles(1);
//
//  } else {
//
//    const int fcsize = _conf.fock_config().size();
//
//    std::set<int> hole_irreps;
//    for (int i = b; i < fcsize; i+=2)
//      if ( _conf[i] == 0 )
//        hole_irreps.insert(_orb_irrep[i/2]);
//
//    for (int i = b; i < fcsize; i+=2)
//      if ( _conf[i] == 1)
//        if ( hole_irreps.count(_orb_irrep[i/2]) )
//          _plist[_np++] = i;
//
//  }

  _np = 0;

  // alpha
  this->count_irreps(_conf.a_unocc_to_site().size(),
                     _nirr,
                     &_orb_irrep[0],
                     &_conf.a_unocc_to_site()[0],
                     &_aa_hole_pair_irrep_count[0]);
  for (std::vector<int>::const_iterator it = _conf.a_occ_to_site().begin(); it != _conf.a_occ_to_site().end(); it++)
    if ( _aa_hole_pair_irrep_count[_orb_irrep[*it]] > 0 )
      _plist[_np++] = 2*(*it)+0;

  // beta
  this->count_irreps(_conf.b_unocc_to_site().size(),
                     _nirr,
                     &_orb_irrep[0],
                     &_conf.b_unocc_to_site()[0],
                     &_bb_hole_pair_irrep_count[0]);
  for (std::vector<int>::const_iterator it = _conf.b_occ_to_site().begin(); it != _conf.b_occ_to_site().end(); it++)
    if ( _bb_hole_pair_irrep_count[_orb_irrep[*it]] > 0 )
      _plist[_np++] = 2*(*it)+1;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  stores all holes matching the supplied particle
///
/// \param[in]      p1       the particle site for which we want matching holes
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::fqmc::DistProposer::get_valid_holes(const int p1) {

  const int p_irrep = _orb_irrep[p1/2];

  const int b = p1 % 2;

  _np = 0;

  if (b == 0) {
    for (std::vector<int>::const_iterator it = _conf.a_unocc_to_site().begin(); it != _conf.a_unocc_to_site().end(); it++)
      if ( _orb_irrep[*it] == p_irrep )
        _plist[_np++] = 2*(*it)+0;
  } else {
    for (std::vector<int>::const_iterator it = _conf.b_unocc_to_site().begin(); it != _conf.b_unocc_to_site().end(); it++)
      if ( _orb_irrep[*it] == p_irrep )
        _plist[_np++] = 2*(*it)+1;
  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  stores all particle pairs with valid matching hole pairs in the pair vector
///
///////////////////////////////////////////////////////////////////////////////////////////////////
void formic::fqmc::DistProposer::get_valid_particle_pairs() {

//  // aa case
//  const int n_aa_holes = _config.a_unocc_to_site().size();
//  if ( n_aa_holes > 1 ) {
//
//    // record all the aa hole pairs' irreps
//    int * irrep_ptr = &_aa_hole_pair_irreps[0];
//    const int * const unocc_ptr = &_config.a_unocc_to_site()[0];
//    for (int i = 1; i < n_aa_holes; i++)
//    for (int j = 0; j < i; j++)
//      *irrep_ptr++ = this->orb_pair_dp(unocc_ptr[i], unocc_ptr[j]);
//
//    // ensure there are not too many hole pairs
//    if ( irrep_ptr - &_aa_hole_pair_irreps[0] > _aa_hole_pair_irreps.size() )
//      throw formic::Exception("too many aa hole pairs!!!");
//
//    // sort the aa hole pair irreps
//    formic::quicksort(&_aa_hole_pair_irreps[0], irrep_ptr, &_quicksort_work[0]);
//
//    // get the number of hole pairs
//    _n_aa_hole_pairs = irrep_ptr - &_aa_hole_pair_irreps[0];
//
//    // place the irreps that are present in the front of the array and count them
//    _n_aa_hole_pair_irrep = 0;
//    for (int i = 1; i < _n_aa_hole_pairs; i++) {
//      // _anti_delta_table[i] = 0  if  i is 0,   1 if i > 0
//      _n_aa_hole_pair_irrep += _anti_delta_table[ _aa_hole_pair_irreps[i] - _aa_hole_pair_irreps[_n_aa_hole_pair_irrep] ];
//      _aa_hole_pair_irreps[_n_aa_hole_pair_irrep] = _aa_hole_pair_irreps[i];
//    }
//    _n_aa_hole_pair_irrep++;
//
//    p = 0;
//    for (int i = 1; i < n; i++)
//      if ( a[i] > a[p] )
//        a[++p] = a[i];
//    p++;
//
//  }

  // get hole pair irreps
  std::set<int> aa_hole_irreps;
  std::set<int> bb_hole_irreps;
  std::set<int> ab_hole_irreps;
  {

    std::set<int> all_irrep;
    for (int i = 0; i < _nirr; i++)
      all_irrep.insert(i);

    // aa hole pair irreps
    _np = 0;
    this->get_pairs(0, 0, 0, all_irrep);
    for (int p = 0; p < _np; p++)
      aa_hole_irreps.insert( this->site_pair_dp(_plist[2*p+0], _plist[2*p+1]) );

    // bb hole pair irreps
    _np = 0;
    this->get_pairs(1, 1, 0, all_irrep);
    for (int p = 0; p < _np; p++)
      bb_hole_irreps.insert( this->site_pair_dp(_plist[2*p+0], _plist[2*p+1]) );

    // ab hole pair irreps
    _np = 0;
    this->get_pairs(0, 1, 0, all_irrep);
    for (int p = 0; p < _np; p++)
      ab_hole_irreps.insert( this->site_pair_dp(_plist[2*p+0], _plist[2*p+1]) );

  }

  // get particle pairs with matching hole pairs
  _np = 0;
  this->get_pairs(0, 0, 1, aa_hole_irreps);
  this->get_pairs(1, 1, 1, bb_hole_irreps);
  this->get_pairs(0, 1, 1, ab_hole_irreps);

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Proposes a configuration difference corresponding to a 1-electron move
///
/// \param[in]      userinp  the user's input options
/// \param[in]      config   the configuration the walker is currently on
/// \param[out]     diff     on exit, the configuration difference object describing the move
/// \param[in,out]  lfg      a random number generator
///
/// \return  the proposal density ratio
///
///////////////////////////////////////////////////////////////////////////////////////////////////
double formic::fqmc::DistProposer::propose_1e(const formic::InputBase & userinp,
                                              const formic::Configuration & config,
                                              formic::ConfigDiff & diff,
                                              formic::LaggedFibonacci<> & lfg) {

  // keep track of how much time is taken in this function
  formic::ScopedTimer stimer(formic::timers::propose_1e.get());
  //formic::start_timer("DistProposer::propose_1e");
  //formic::timers::propose_1e.start();

  // store the configuration
  _conf = config;

  // get the particles that can be moved
  this->get_valid_particles();

  // if there are no possible single moves, perform a double move instead
  if ( _np == 0 ) {
    //formic::stop_timer("DistProposer::propose_1e");
    //formic::timers::propose_1e.stop();
    stimer.clear(); // stops the timer before running 2e proposal function
    //return this->propose_2e(userinp, config, diff, lfg);
    return this->fast_propose_2e(userinp, config, diff, lfg);

    //std::cout << "offending configuration:    ";
    //for (int i = 0; i < _conf.fock_config().size(); i+=2)
    //  std::cout << _conf[i];
    //std::cout << "  ";
    //for (int i = 1; i < _conf.fock_config().size(); i+=2)
    //  std::cout << _conf[i];
    //std::cout << "      _orb_irrep:";
    //for (int i = 0; i < _orb_irrep.size(); i++)
    //  std::cout << _orb_irrep[i];
    //std::cout << std::endl;
    //throw formic::Exception("DistProposer::propose_1e found no valid forward particles");
  }

  // remember the number of possible forward particles
  const int nfp = _np;

  // choose a particle at random
  const int p1 = _plist[lfg() % _np];

  // get the valid holes
  this->get_valid_holes(p1);
  if ( _np == 0 )
    throw formic::Exception("DistProposer::propose_1e found no valid forward holes");

  // get the distance weights
  this->fast_get_distance_weights(p1);

  // pick a hole at random
  int cp = -1;
  for (double rn = formic::uni_01(lfg); rn >= 0.0 && cp+1 < _np; )
    rn -= _dw[++cp];
  if ( cp < 0 || cp >= _np )
    throw formic::Exception("DistProposer::propose_1e has bad index for forward hole");

  // get the hole
  const int h1 = _plist.at(cp);

  // get the forward proposal probability
  const double fpp = _dw.at(cp) / nfp;

  //if ( fpp == 0.0 )
  //  throw formic::Exception("DistProposer::propose_1e has chosen a zero forward proposal probability");

  // create the difference object
  diff.newly_full().resize(1);
  diff.newly_full()[0] = h1;
  diff.newly_empty().resize(1);
  diff.newly_empty()[0] = p1;

  // if a zero probability proposition was chosen, return zero to ensure the move is rejected
  if ( fpp == 0.0 )
    return 0.0;

  // update the configuration object according to the chosen move
  _conf.update(diff);

  // prepare the array of possible reverse move particles
  this->get_valid_particles();
  if ( _np == 0 )
    throw formic::Exception("DistProposer::propose_1e found no valid reverse particles");

  // remember the number of reverse particles
  const int nrp = _np;

  // ensure the reverse particle is present
  for (cp = 0; cp < _np; cp++) {
    if ( _plist[cp] == h1)
      break;
  }
  if ( cp >= _np )
    throw formic::Exception("DistProposer::propose_1e cannot find forward hole in reverse particles.");

  // get the valid reverse holes
  this->get_valid_holes(h1);
  if ( _np == 0 )
    throw formic::Exception("DistProposer::propose_1e found no valid reverse holes");

  // get the reverse distance weights
  this->fast_get_distance_weights(h1);

  // find the reverse move and return the proposal density ratio
  for (cp = 0; cp < _np; cp++) {
    if (_plist[cp] == p1) {
      //formic::stop_timer("DistProposer::propose_1e");
      //formic::timers::propose_1e.stop();
      const double rpp = _dw.at(cp) / nrp; // reverse proposal probability
      return ( rpp / fpp ); // return the ratio of reverse to forward proposal probabilities
    }
  }

  //formic::stop_timer("DistProposer::propose_1e");
  //formic::timers::propose_1e.stop();

  throw formic::Exception("DistProposer::propose_1e failed to find reverse move");

  return 1.0;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Proposes a configuration difference corresponding to a 2-electron move
///
/// \param[in]      userinp  the user's input options
/// \param[in]      config   the configuration the walker is currently on
/// \param[out]     diff     on exit, the configuration difference object describing the move
/// \param[in,out]  lfg      a random number generator
///
/// \return  the proposal density ratio
///
///////////////////////////////////////////////////////////////////////////////////////////////////
double formic::fqmc::DistProposer::fast_propose_2e(const formic::InputBase & userinp,
                                                   const formic::Configuration & config,
                                                   formic::ConfigDiff & diff,
                                                   formic::LaggedFibonacci<> & lfg) {

  // keep track of how much time is taken in this function
  formic::ScopedTimer stimer(formic::timers::propose_2e.get());
  //formic::start_timer("DistProposer::fast_propose_2e");
  //formic::timers::propose_2e.start();

  // store the configuration
  _conf = config;

  // determine all the hole and particle pairs and their irreps for the forward move
  this->prep_pairs_and_irreps();

  // based on the hole pairs' irreps, place the allowed particle pairs in _plist
  _np = 0;
  _np += this->get_irrep_allowed_pairs(0, 0, _n_all_aa_part_pairs, _nirr, &_all_aa_part_pairs[0], &_aa_hole_pair_irrep_count[0], &_plist[0] + 3*_np);
  _np += this->get_irrep_allowed_pairs(1, 1, _n_all_bb_part_pairs, _nirr, &_all_bb_part_pairs[0], &_bb_hole_pair_irrep_count[0], &_plist[0] + 3*_np);
  _np += this->get_irrep_allowed_pairs(0, 1, _n_all_ab_part_pairs, _nirr, &_all_ab_part_pairs[0], &_ab_hole_pair_irrep_count[0], &_plist[0] + 3*_np);
  if ( _np == 0 )
    throw formic::Exception("DistProposer::fast_propose_2e found no valid forward particle pairs");

  // remember the number of possible forward proposal pairs
  const int nfp = _np;

  // choose a pair of occupied orbitals at random
  const int forward_pair = lfg() % _np;
  const int fwd_irr = _plist.at(3*forward_pair+0);
  const int p1 = _plist.at(3*forward_pair+1);
  const int p2 = _plist.at(3*forward_pair+2);
  std::fill(_chosen_pair_irrep_counts.begin(), _chosen_pair_irrep_counts.end(), 0);
  _chosen_pair_irrep_counts.at(fwd_irr) = 1;

  // get the valid forward move hole pairs
  if      ( p1 % 2 == 0 && p2 % 2 == 0 )
    _np = this->get_irrep_allowed_pairs(0, 0, _n_all_aa_hole_pairs, _nirr, &_all_aa_hole_pairs[0], &_chosen_pair_irrep_counts[0], &_plist[0]);
  else if ( p1 % 2 == 1 && p2 % 2 == 1 )
    _np = this->get_irrep_allowed_pairs(1, 1, _n_all_bb_hole_pairs, _nirr, &_all_bb_hole_pairs[0], &_chosen_pair_irrep_counts[0], &_plist[0]);
  else if ( p1 % 2 == 0 && p2 % 2 == 1 )
    _np = this->get_irrep_allowed_pairs(0, 1, _n_all_ab_hole_pairs, _nirr, &_all_ab_hole_pairs[0], &_chosen_pair_irrep_counts[0], &_plist[0]);
  else
    throw formic::Exception("DistProposer::fast_propose_2e did not expect a forward hole pair in ba order");
  if ( _np == 0 )
    throw formic::Exception("DistProposer::fast_propose_2e found no valid forward hole pairs");

  // get the distance weights
  this->fast_get_distance_weights(p1, p2);

  // pick a hole pair at random
  int cp = -1; // chosen pair
  for (double rn = formic::uni_01(lfg); rn >= 0.0 && cp+1 < _np; )
    rn -= _dw[++cp];
  if ( cp < 0 || cp >= _np )
    throw formic::Exception("DistProposer::fast_propose_2e has bad index for forward hole pair");

  // get the holes
  const int h1 = _plist.at(3*cp+1);
  const int h2 = _plist.at(3*cp+2);

  // get the forward proposal probability
  const double fpp = _dw.at(cp) / nfp;

  //if ( fpp == 0.0 )
  //  throw formic::Exception("DistProposer::fast_propose_2e has chosen a zero forward proposal probability");

  // create the difference object
  diff.newly_full().resize(2);
  diff.newly_full()[0] = h1;
  diff.newly_full()[1] = h2;
  diff.newly_empty().resize(2);
  diff.newly_empty()[0] = p1;
  diff.newly_empty()[1] = p2;

  // if a zero probability proposition was chosen, return zero to ensure the move is rejected
  if ( fpp == 0.0 ) {
    //formic::timers::propose_2e.stop();
    return 0.0;
  }

  // update the configuration object according to the chosen move
  _conf.update(diff);

  // determine all the hole and particle pairs and their irreps for the reverse move
  this->prep_pairs_and_irreps();

  // based on the hole pairs' irreps, place the allowed particle pairs in _plist
  _np = 0;
  _np += this->get_irrep_allowed_pairs(0, 0, _n_all_aa_part_pairs, _nirr, &_all_aa_part_pairs[0], &_aa_hole_pair_irrep_count[0], &_plist[0] + 3*_np);
  _np += this->get_irrep_allowed_pairs(1, 1, _n_all_bb_part_pairs, _nirr, &_all_bb_part_pairs[0], &_bb_hole_pair_irrep_count[0], &_plist[0] + 3*_np);
  _np += this->get_irrep_allowed_pairs(0, 1, _n_all_ab_part_pairs, _nirr, &_all_ab_part_pairs[0], &_ab_hole_pair_irrep_count[0], &_plist[0] + 3*_np);
  if ( _np == 0 )
    throw formic::Exception("DistProposer::fast_propose_2e found no valid reverse particle pairs");

  // remember the number of reverse particle pairs
  const int nrp = _np;

  // get the valid reverse move hole pairs
  if      ( p1 % 2 == 0 && p2 % 2 == 0 )
    _np = this->get_irrep_allowed_pairs(0, 0, _n_all_aa_hole_pairs, _nirr, &_all_aa_hole_pairs[0], &_chosen_pair_irrep_counts[0], &_plist[0]);
  else if ( p1 % 2 == 1 && p2 % 2 == 1 )
    _np = this->get_irrep_allowed_pairs(1, 1, _n_all_bb_hole_pairs, _nirr, &_all_bb_hole_pairs[0], &_chosen_pair_irrep_counts[0], &_plist[0]);
  else if ( p1 % 2 == 0 && p2 % 2 == 1 )
    _np = this->get_irrep_allowed_pairs(0, 1, _n_all_ab_hole_pairs, _nirr, &_all_ab_hole_pairs[0], &_chosen_pair_irrep_counts[0], &_plist[0]);
  else
    throw formic::Exception("DistProposer::fast_propose_2e did not expect a reverse hole pair in ba order");
  if ( _np == 0 )
    throw formic::Exception("DistProposer::fast_propose_2e found no valid reverse hole pairs");

  // get the reverse move distance weights
  this->fast_get_distance_weights(h1, h2);

  // find the reverse move and get the proposal density ratio
  //std::cout << boost::format(" forward move:\n%4i %4i %4i") % fwd_irr % p1 % p2 << std::endl;
  for (int x = 1, three_np = 3*_np; x < three_np; x += 3) {
    const bool found =    ( _plist[x] == p1 && _plist[x+1] == p2 )
                       || ( _plist[x] == p2 && _plist[x+1] == p1 );
    if (found) {
      //formic::stop_timer("DistProposer::fast_propose_2e");
      //formic::timers::propose_2e.stop();
      const double rpp = _dw.at(x/3) / nrp; // reverse proposal probability
      return ( rpp / fpp ); // return the ratio of reverse to forward proposal probabilities
    }
  }

  //formic::stop_timer("DistProposer::fast_propose_2e");
  //formic::timers::propose_2e.stop();

  throw formic::Exception("DistProposer::fast_propose_2e failed to find reverse move");

  return 1.0;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Proposes a configuration difference corresponding to a 2-electron move
///
/// \param[in]      userinp  the user's input options
/// \param[in]      config   the configuration the walker is currently on
/// \param[out]     diff     on exit, the configuration difference object describing the move
/// \param[in,out]  lfg      a random number generator
///
/// \return  the proposal density ratio
///
///////////////////////////////////////////////////////////////////////////////////////////////////
double formic::fqmc::DistProposer::propose_2e(const formic::InputBase & userinp,
                                              const formic::Configuration & config,
                                              formic::ConfigDiff & diff,
                                              formic::LaggedFibonacci<> & lfg) {

  // keep track of how much time is taken in this function
  formic::ScopedTimer stimer(formic::timers::propose_2e.get());
  //formic::start_timer("DistProposer::propose_2e");

  // store the configuration
  _conf = config;

  // get the valid particle pairs
  this->get_valid_particle_pairs();

  if ( _np == 0 )
    throw formic::Exception("DistProposer::propose_2e found no valid forward particle pairs");

  // remember the number of possible forward proposal pairs
  const int nfp = _np;

  // choose a pair of particles at random
  const int forward_pair = lfg() % _np;
  const int p1 = _plist[2*forward_pair+0];
  const int p2 = _plist[2*forward_pair+1];

  // get the valid hole pairs
  {
    std::set<int> irrep_set;
    irrep_set.insert( this->site_pair_dp(p1, p2) );
    _np = 0;
    this->get_pairs(p1 % 2, p2 % 2, 0, irrep_set);
  }

  if ( _np == 0 )
    throw formic::Exception("DistProposer::propose_2e found no valid forward hole pairs");

  // get the distance weights
  this->get_distance_weights(p1/2, p2/2);

  // pick a hole pair at random
  int cp = -1; // chosen pair
  for (double rn = formic::uni_01(lfg); rn >= 0.0 && cp+1 < _np; )
    rn -= _dw[++cp];

  if ( cp < 0 || cp >= _np )
    throw formic::Exception("DistProposer::propose_2e has bad index for forward hole pair");

  // get the holes
  const int h1 = _plist.at(2*cp+0);
  const int h2 = _plist.at(2*cp+1);

  // get the forward proposal probability
  const double fpp = _dw.at(cp) / nfp;

  //if ( fpp == 0.0 )
  //  throw formic::Exception("DistProposer::propose_2e has chosen a zero forward proposal probability");

  // create the difference object
  diff.newly_full().resize(2);
  diff.newly_full()[0] = h1;
  diff.newly_full()[1] = h2;
  diff.newly_empty().resize(2);
  diff.newly_empty()[0] = p1;
  diff.newly_empty()[1] = p2;

  // if a zero probability proposition was chosen, return zero to ensure the move is rejected
  if ( fpp == 0.0 )
    return 0.0;

  // update the configuration object according to the chosen move
  _conf.update(diff);

  // prepare the array of possible reverse move particle pairs
  this->get_valid_particle_pairs();

  if ( _np == 0 )
    throw formic::Exception("DistProposer::propose_2e found no valid reverse particle pairs");

  // remember the number of reverse particle pairs
  const int nrp = _np;

  // ensure the reverse particle pair is present
  for (cp = 0; cp < _np; cp++) {
    bool found =    ( _plist[2*cp+0] == h1 && _plist[2*cp+1] == h2 )
                 || ( _plist[2*cp+0] == h2 && _plist[2*cp+1] == h1 );
    if (found)
      break;
  }
  if ( cp >= _np )
    throw formic::Exception("DistProposer::propose_2e cannot find forward holes in reverse particles.");

  // get the valid reverse hole pairs
  {
    std::set<int> irrep_set;
    irrep_set.insert( this->site_pair_dp(h1, h2) );
    _np = 0;
    this->get_pairs(h1 % 2, h2 % 2, 0, irrep_set);
  }

  if ( _np == 0 )
    throw formic::Exception("DistProposer::propose_2e found no valid reverse hole pairs");

  // get the reverse distance weights
  this->get_distance_weights(h1/2, h2/2);

  // find the reverse move and get the proposal density ratio
  for (cp = 0; cp < _np; cp++) {
    bool found =    ( _plist[2*cp+0] == p1 && _plist[2*cp+1] == p2 )
                 || ( _plist[2*cp+0] == p2 && _plist[2*cp+1] == p1 );
    if (found) {
      //formic::stop_timer("DistProposer::propose_2e");
      const double rpp = _dw.at(cp) / nrp; // reverse proposal probability
      return ( rpp / fpp ); // return the ratio of reverse to forward proposal probabilities
    }
  }

  //formic::stop_timer("DistProposer::propose_2e");

  throw formic::Exception("DistProposer::propose_2e failed to find reverse move");

  return 1.0;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief  Proposes a configuration difference corresponding to a possible walker move
///
/// \param[in]      userinp  the user's input options
/// \param[in]      config   the configuration the walker is currently on
/// \param[out]     diff     on exit, the configuration difference object describing the move
/// \param[in,out]  lfg      a random number generator
///
/// \return  the proposal density ratio
///
///////////////////////////////////////////////////////////////////////////////////////////////////
double formic::fqmc::DistProposer::propose(const formic::InputBase & userinp,
                                           const formic::Configuration & config,
                                           formic::ConfigDiff & diff,
                                           formic::LaggedFibonacci<> & lfg) {

  const double rn = formic::uni_01(lfg);
  if ( _2e_prob == 0.0 )
    return this->propose_1e(userinp, config, diff, lfg);
  if ( rn < _2e_prob || _2e_prob == 1.0 )
    //return this->propose_2e(userinp, config, diff, lfg);
    return this->fast_propose_2e(userinp, config, diff, lfg);
  return this->propose_1e(userinp, config, diff, lfg);

}
